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ABSTRACT 

We present a finite difference code intended for computing linear, adiabatic, 
nonradial pulsations of spherical stars. This code is based on a slight 
modification of the general Newton - Raphson technique in order to handle the 
relaxation of the eigenvalue (square of the eigenfrequency) of the modes and 
their corresponding eigenf unctions. This code has been tested computing the 
pulsation spectra of polytropic spheres finding a good agreement with previous 
work. 

Then, we have coupled this code to our evolutionary code and applied it to 
the computation of the pulsation spectrum of a low mass, pure - helium white 
dwarf of 0.3 Mq for a wide range of effective temperatures. In making this 
calculation we have taken an evolutionary time step short enough such that 
eigenmodes corresponding to a given model are used as initial approximation to 
those of the next one. Specifically, we have computed periods, period spacing, 
eigenfunctions, weight functions, kinetic energies and variational periods for 
a wide range of modes. To our notice this is the first effort in studying the 
pulsation properties of helium white dwarfs. The solution we have found 
working with these realistic white dwarf models are in good accord with the 
predictions of the asymptotic theory of Tassoul (1980) for high order modes. 
This indicates that the code presented here is able to work adequately also with 
reahstic stellar models. 

Subject headings: Methods: numerical - stars: Interiors - stars: oscillations - 
stars: white dwarfs 
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1. INTRODUCTION 

At present, nonradial stellar pulsations constitute a valuable tool in extracting 
information about the properties and internal structure of objects that undergo such kind 
of oscillations. As observational techniques are refined, a growing variety of stellar objects 
is discovered that undergo nonradial oscillations. The most extensively studied case is our 
Sun, which has been shown to vibrate in a large amount of modes, investigation of which has 
been so powerful in revealing clues of its internal structure that it is known as the discipline 
of Helioseismology. Also, it has been established that different kinds of stellar objects which 
are located in a variety of places in the HR diagram, covering several evolutionary stages, 
undergo nonradial pulsations, e.g.: Ap, WR, 5 Scuti, jS Cephei, and variable white dwarfs 
(hereafter WDs): DAV, DBV, DOV and PNNV. The study of oscillations in such variable 
stars, and the subsequent comparison and fitting with theoretical pulsation patterns is now 
known as A stero seismology. For a general introduction to this topic see, e.g.. Cox (1980); 
Unno, et al. (1989); Brown & Gilliland (1994); and Gaustchy & Saio (1995; 1996). 

Prom the observational point of view (with the obvious exception of our Sun) , WDs 
represent the best established kind of nonradial pulsators. WDs are the final fate for low 
and intermediate mass stars. During its cooling, they go across very narrow instability 
strips inside which suffer nonradial pulsations in the branch of g- (gravity) modes with 
periods between 100 and 2000 sec and amphtudes up to 0.3 magnitudes. 

As WDs cools down, the outermost layers main constituent (hydrogen in DAVs, 
helium in DBVs, and probably carbon or oxygen in DOVs and PNNVs; see e.g. Brown & 
Gilliland 1994 and Gaustchy & Saio 1995; 1996) gets partial ionization conditions which 
forces the development of an outer convection zone (OCZ). In the bottom of such zone the 
instabilities are originated {k mechanism of overstability) . This gives rise to the blue edge 
of the instability strip (see, e.g., Dziembowski & Koester 1981 and Winget, et al. 1982; 



-4- 



see also Tassoul, Fontaine & Winget 1990, hereafter TFW)[|. The physical reason for the 
appearance of a red edge is still under debate. 

At present, there is a general consensus that variable WDs are very interesting objects 
for asteroseismological studies. Their very simple internal structure allows us to predict 
theoretically the modes with a very high degree of detail and sophistication. Also, they 
have a very rich spectrum of periods which may give us information about stellar masses 
(by means of measuring the mean period spacing), core composition (by measuring the rate 
of change of the period P), mass of the surface helium and hydrogen layers (if present) 
(by means of the departures from uniform period spacing and mode trapping), velocity of 
rotation and strength of the magnetic field (studying the fine structure of a multiplet), 
etc. (see, e.g., Bradley 1996). So, it is not surprising that in the last years DAVs and 
DBVs as well as DOVs have been the preferred target for the network called "Whole Earth 
Telescope" (WET). WET observations have been of an unprecedent quality which in some 
cases (see, e.g., Winget et al. 1991; Kawaler 1993; Bradley & Winget 1994) allowed to 
apply the powerful tools above mentioned. 

In view of the growing importance of this line of research, in our Observatory we 
have begun the study of nonradial oscillations. To this end, we have developed a code to 
compute eigenfunctions and eigenfrequencies for spherical objects in the linear, adiabatic 
approximation. The main purpose of the present paper is to describe the numerical 
techniques employed in our program to search for and to compute the modes and its 
application to the case of a 0.3 Mq pure helium WD. This code is based on a modification 
of the general Newton - Raphson technique presented in Kippenhahn, Weigert & Hofmeister 
(1967) to solve the set of difference equations which represent the differential equations of 

^Recently Goldreich & Wu (1999) have proposed the "convective driving mechanism" as 
the responsible for overstability of g-modes in DA WDs. 
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linear, adiabatic, nonradial pulsations of spherical stars. 

In order to test the code we have applied it to the case of polytropic spheres whose 
eigenfrequencies have been computed with very high accuracy (Christensen - Dalsgaard & 
Mullan 1994). Also, we have coupled our pulsation code to an stellar evolution code. We 
have done it in order to be able to follow the changes in the pulsational spectrum of a given 
model during its evolution with a number of steps as large as it would be needed by specific 
studies. In making such computation we selected an evolutionary time step short enough 
such that eigenvalues and eigenf unctions corresponding to a given model are used as a 
initial approximation to those of the next one. This is very useful for a relaxation code like 
this and is intended to work as it is the case of the Heneyey technique in stellar evolution. 

The rest of the paper is organized as follows: Section |^ is devoted to present the 
differential equations of linear, adiabatic, nonradial pulsations we have to solve. In Section ^ 
a detailed description of the general Newton - Raphson technique we employ for computing 
the eigenmodes is presented. Section ^ presents the results for the case of polytropic 
spheres, and in Section ^ we describe the work we have done with a realistic helium WD 
model. Finally, in Section ^ summarize our results. 

2. The Differential Equations of Nonradial Pulsations 

The differential equations that govern linear, nonradial pulsations of spherical stars in 
the adiabatic approximation are (Unno, et al. 1989): 
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X 



-£ - (C, a;2 - A*) y, + (A*-U+l)y,- A* y^, 



(2) 



X 



dy 



(3) 



x^^U A*y, + UVgy2 + 



e{e + i)-uVg 



yz-U y4. 



(4) 



X is the independent variable defined as x — r / R where r is the radial coordinate and R is 

the stellar radius. The dependent variables are defined as 



yi = — , 

r 



(5) 



(6) 



and 



9 r 



(7) 



1 d<^ 



y4 = - 



dr 



(8) 



a; 



(9) 



represents the radial displacement of the fluid and p and $' are the Eulerian variations 
of the equilibrium values of the pressure p and the gravitational potential $ respectively. 



- 7- 



g is the local acceleration of gravity and is the stellar mass, uj"^ is the dimensionless 
square of the angular frequency of oscillation, cr^. 

Among the coefficients of Eqs. there appears the harmonic degree corresponding 
to the spherical harmonic „j(^, 0) which describes the angular dependence of the oscillation 
pattern. The adimensional quantities V^, f/, Ci, and A* , inherent to the non - perturbed 
model, are defined as: 



_Y_ _ _ ]_ ^^^P _ ^ _ 9r p , . 

' Ti Ti dlnr Ti p ' ^ ^ 



dlnr M. ' ^ ^ 



r,^n /l(ilnp(ilnp\ 
A* = -N^ = r[-—^-—^]. (13) 

Here Fi is the first adiabatic index, is the square of the local velocity of sound, is the 
mass contained in a (non - perturbed) sphere of radius r, and N"^ is the Brunt - Vaisala 
frequency. 



Eqs. (|T1-|D together with the adequate boundary conditions conform the eigenvalue 
problem to be solved. 
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3. The Numerical Code 

In order to write the pulsation equations in a form adequate for numerical calculations 
we divide our non - perturbed model in a number of spherical, concentric shells. This is 
equivalent to divide the dominion of the independent variable in mesh points (A^ — 1 
shells) non - necessarily equidistantly spaced Xj,j = 1, ■ ■ ■ , N. In our treatment we define 
xi = 1 as the surface and xn = as the central point of the model. 

Now, we replace Eqs. (0-0) by difference equations. The system of equations may be 
rewritten as 

^ = /i(l/i,l/2,l/3,2/4, A), i = l,2,3,4, (14) 

where A = a;^. In finite differences, Eq. (|I3|) may be written asQ 



[yi]j+i - [yi]j ^ 

Xj+i — Xj 

where, for any quantity ip, 



ft ([2/i,2/2,2/3,2/4]j+i; a) ;i = l,2,3,4;j = l,---,iV- 1 (16) 



= o ' (1^) 



^We have also employed finite differences of the form 



_ = I (/. ([2/1,1/2, 2/3, 2/4],; A) + /, ([yi, 7/2, 2/3, 2/4],+i; A)) ; z = 1, 2, 3, 4; J = 1, ■ ■ ■ , N-1. 

(15) 

Applying such form for the difference equations we have found results negligibly different 
from those presented in this paper. 
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being the value of the quantity ip at the meshpoint j. The outer boundary conditions 
are given byQ (Unno, et ah 1989) 

bi]i - [y2\i + [1/3] 1 = 0, 

(£ + 1) [ysli + [2/4]i = 0, (18) 
[yi]i = 1. 

The last equation is the normalization condition, usually employed in previous works. The 
central boundary conditions are given by (Unno, et al. 1989) 

[l/iU ^ - b]^ = 0, 

(i9j 

^ [VsIn - bJ^v = 0. 

In order to solve this system of difference equations we shall use the general Newton - 
Raphson technique closely following the formulation presented in Kippenhahn, et al. (1967) 
for the case of stellar evolution. We shall assume we are able to get an approximate solution 
for the system and we want to improve it iteratively. In the case that this initial solution 
is not far from the actual solution of the system, we may expand the equations up to first 
order in the corrections for the values of the eigenfunctions at each meshpoint and also for 
the eigenvalue (square of the eigenfrequency) . In this standard way we get a linear system 
of equations in the corrections that must be solved. 

In condensed way, the algebraic system of equations for the first order corrections may 
be expressed as 

^The first two of these equations are the so called "zero boundary condition" , which we 
shall employ in the treatment of polytropes. In the case of realistic WD models we have 
replaced the first one by Eq. (|30|) (see below). 
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5[y,V + ■ ■ ■ + 1^ ^Mi + ^ = -B,- A; = 1, 2, 3, (20) 



d[yi]i d[y4\i dX 



^ 5A = -Gi; ^ = 1,2,3,4; j = 1, 2, . . . , iV - 1, 



(21) 



^^"^ 5bi]iv + ■ ■ ■ + 1% 5[?/4]^ + ^ = -C^; m = 1, 2, (22) 



d[yi]N d[yA]N dX 

where S[yi]j are small corrections to the eigenfunction yi at the j— mesh point, and 6X 
stand for the correction to the eigenvalue A. Bk, G\, and Cm are the values of the difference 
equations when applied the solution to be iteratively improved. Obviously, all of them must 
be zero when evaluated with their exact solution. 

Now, we have to invert a big matrix, that has non - zero elements only in blocks 
located on the diagonal and in the last column (because of the derivative with respect to 
the eigenvalue). Notice that this is an important difference compared with the case of 
stellar evolution. Thus, in handling the big matrix we need a specific algorithm to solve for 
the corrections. 



Evaluating Eq. (20) for k = 1, 2, 3 and Eq. (|2l|) for i = 1, 2, 3, 4 and j = 1, we may 



write the first block of such matrix as 
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dBi 
dlyi]i 


dBi 
d[y4]i 




















dBi 
d\ 


-B, 






dB2 


dB2 

d[y4.]i 













S[y2\i 







dB2 
d\ 


-B2 






dB:i 
9[2/i]i 


dB3 
d[y4]i 













S[y3]i 







dBa 
dX 


-Bs 






dG] 
d[yi]i 


dG] 
d[y4]i 


dG] 
d[yi]2 


dG] 
d[y2]2 


dG] 
dlysh 




S[y4\i 
S[yih 

(5 [1/2] 2 




dG] 

d[y4]2 


dG] 
dX 


-G\ 




5\ 
1 


OGl 
dlyi]i 


OGl 
d[y4.]i 


dGl 
d[yi]2 


dGl 
d[y2]2 


dGl 
dlysh 




5 [1/3] 2 




dGl 
d[y4]2 


dGl 
d\ 


-G\_ 







(23) 

Now, we define vectors Ug, Vg, Wg, {s — 1, - ■ ■ , AN — 5) such that we may write 



S[yi]i 




S[y2]i 




S[y3]i 




S[y4\i 




^[yih 




S[y2]2 




Siysh 





U2 V2 W2 



U7 V7 W7 



S[y4h 

6X 
1 



(24) 



for this block. The coefficients of this vectors may be calculated easily solving 
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dlyi]i 


d[yA\i 











832 
dlyi]i 


dB2 
d[y4\i 











dB3 
dlyi]i 


dBa 
d[y4]i 











dG] 


dG] 


dG] 


dG] 


dG] 


dlyi]i 


dly4]i 


d[yi]2 


d[y2]2 


dlysh 


dGl 


dGl 


dGl 


dGl 


dGl 


dlyi]i 


dly4]i 


d[yi]2 


d[y2]2 


dlysh 



Ui Vi Wi 

U2 V2 W2 



U7 V7 W: 
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dG] 

dlyih 



dGl 

d[y4,]2 



dBi 
' d\ 

dB2 
d\ 

dBi 
d\ 

dG] 
' dX 



dGl 
' dX 



-B2 



(25) 



From Eq. (^Tj), and after short manipulation, we may write an arbitrary block of the big 
matrix (except the central) in the form 



where j = 2, 
the system 







f/4, 






















Vij+2 








U4j+3 







6X 
1 



(26) 



N — 2. Now, the coefficients of vectors Us,Vs, Wg may be evaluated solving 



"3 



dG{ 



d[y 



dG] 



dGi 



d[y: 



'2\j+l 



d[yi]j+i d[y: 

d[yi]j+i 
dGj 



dGi 



d[yi]j+i d[y2]j 



'2\j+l 

d[y2]j+i 
dGj 



dG] 
d[ya]j+i 
dGj 

d[yi\3+i 

dGj 
d[y3]j+i 















UAj+2 


V4j+2 


W4j+2 






W^4i+3 



dG{ 

d[y4,h+i 

dGj 
d[y4]j+i 

dGj 
d[y4]j+i 

dGl 
d[y4]j+i 



-Pi 
-Pi 
-Pi 
-Pi 



-li 
-73 
-74 



(27) 



where quantities a\^(5l^^l 



are defined as 
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dGi jj dGi jj dGl jj dG\ 



9G' 



(28) 



with i = 1,2,3,4. 

Finally, in the central point of the model, we have j = — 1 in Eq. (^) and m = 1, 2 
in (P^l). It is easy to show that the last block of the big matrix may be written as 



a 



N-l 



a 



N-l 
3 

at 





1 




9Gf-i 


9Gf-i 


9Gf-i 




d[yi]N 


9[s/2]jv 




'9[y4]jv 


1 




9G^-i 


9G^-i 


9G^-i 




d[yi]N 


9[s/2]jv 




<9[j/4]jV 


1 


9G,f-i 


9G^-i 


9G.^-i 


9G^-i 




d[yi]N 


9[s/2]jv 




<9[y4]jv 


1 


9Gf-i 


aGf-i 


9Gf-i 


9Gf-i 




d[yi]N 


9[s/2]jv 


9[s/3]]v 


d[y4,]N 




dCi 


9Gi 


9Gi 


dCi 




d[yi]N 


9[s/2]jv 


d[y'i]N 


d[y4]N 




dC2 


9G2 


dC2 


dC2 




d[yi\N 


9[s/2]jv 


d[y3]N 


d[y4]N 











-it' 






S[yi]N 




-it' 






S[y2]N 




-it' 






^[yslN 




-it' 


aci 

d\ 




S[y4]N 






dC2 
dX 




6X 




-C2 



In this case the quantities ^, pj^ 'yf ^ are evaluated from Eqs. 



(29) 



at j = - 1. 



We note that Eq. (p9|) can be solved in order to obtain the corrections for eigenfunctions 
yi, 1/2, 2/3 and 2/4 in the central point of the object, and for the eigenvalue A. Also, it must 
be noted that the correction in eigenfunction 7/4 pertaining to the immediate outer mesh 
point of grid is obtained. This correction, namely S[yi]j\f_i, serve as "coupling" between 
the points A^ and A^ — 1. In fact, we may use Eq. (^) with j = A^ — 2 in order to obtain 
the rest of corrections belonging to the eigenfunctions at the mesh point A^ — 1. Next, 
the employment of this procedure for successive downwards values of j (in Eq. |2^) using 
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5[y4]j+i as coupling between the quantities belonging to consecutive meshpoints, as well as 
Eq. (p^, leads to find the corrections for the eigenvalue and eigenf unctions for the complete 
model. These corrections are applied to the initial solution and the algorithm is employed 
iteratively up to the stage at which all the (absolute value) of the corrections are smaller 
than some previously adopted value. At this point, the complete set of difference equations 
that represent Eqs. has been solved. 

In order to search for the first approximation to the eigenfunctions and the eigenvalue 
of a mode, we have applied the method of the discriminant presented in Unno et al. (1989). 
Specifically, the expression adopted reads D{lj) = (i + l) [ysli + [vaIi- Notice that D{ij) = 
corresponds exactly to the second outer boundary condition given in Eqs. ([T^). We refer 
the reader to that book for more details. 

Obviously, for applying this technique we need to define a grid of meshpoints. In 
getting an appropriate distribution of meshpoints we have employed a simple recipe. As 
a first approximation to the solution of the oscillation equations we have computed the 
eigenfunctions and eigenvalue taking the points at which the equilibrium model is defined. 
After convergence, we have redistributed the grid asking for the same criterium usually 
employed in stellar evolution (Kippenhahn, et al. 1967): that the relative variation of 
each eigenfunction inside a zone must be below some prescribed limit (if the eigenfunctions 
are below some other very small, prescribed value we have asked for absolute differences 
in place of relative ones). If necessary, our program add or remove meshpoints, using 
spline interpolation over quantities of non - perturbed model, and lineal interpolation in 
eigenfunctions, since these will be improved after in the general Newton - Raphson stage. 
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4. Application to Polytropic Spheres 

As a first test for our code, we liave applied it to the well - known case of nonradial 
pulsations of polytropic spheres. Polytropes are particularly adequate for the purpose of 
testing, because in this case we separate the inaccuracies of the non - perturbed model 
from those inherent to the nonradial pulsation code. In particular we have computed 
the p- (pressure), f- (fundamental) and g- (gravity) modes of polytropes with indices 
n = 1.5, 2, 2.5, 3, 3.5 and 4 for the case of £ = 2, 3 and 4. 

Because in this case Vg and A* are divergent at surface (see Eqs. |10| and |13]), we have 
integrated the structure by means of an accurate Runge - Kutta technique (Press, et al. 
1992) and most of the meshpoints were located near surface. Here we assumed Fi = |. 

To compare our eigenvalues with those available in the literature we shall employ the 
work of Christensen - Dalsgaard & MuUan (1994) which presents accurate eigenfrequency 
tabulations for a variety of polytropic spheres. For such purpose we have plotted in Fig. 
|1| the relative difference Icu^ — ^ch~D\/^ch~D ^^r the cases of p-modes for n = 1.5,3 and 
i = 2,3, g-modes for n = 3 and i = 2,3, 4, and also for p-modes for n = 4 and i = 2,3. 
The comparison indicates that the higher radial order the larger the relative differences 
are, without any dependency on i value. The absolute value of relative differences found in 
the eigenvalues are ^ 4 x 10~^, except in the n = 3 g-modes, for which the difference is 
^ 2 X 10"'^. As Christensen - Dalsgaard & Mullan (1994) state that their computations 
are accurate to ~ 10~^ we conclude that our code is able to compute the eigenvalues of 
polytropes with a precision of ~ 10^'^ which is enough for our purposes. 
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5. Application to Realistic Models: Nonradial g-Modes in Helium White 

Dweirfs 

With the aim of investigating the behaviour of the code when it is apphed to reahstic 
models, in this section we shall consider nonradial g-modes in evolutionary models of low 
mass, pure hehum WDs. This election is due mainly because of two reasons. First, at 
present we have available a detailed and updated code we have employed in the calculation 
of the evolution of such stars (see, e.g., Althaus & Benvenuto 1997; Benvenuto & Althaus 
1998). Such evolutionary code is fully described in the above-cited papers. Second, the 
computation of nonradial modes in pure helium WDs provide us an important way for 
testing the pulsation code when applied to realistic models. At present, to our notice 
there is no study available in the literature on the pulsation properties of these objects. 
In evaluating our results we shall pay special attention to the asymptotic behaviour of 
periods of oscillation, comparing in particular the spacing of adjacent periods (for the same 
£) with that predicted by asymptotic theory of Tassoul (1980) in stars with homogeneous 
composition. This theory has been employed as check for the nonradial spectrum of 
polytropes (MuUan & Ulrich 1988; MuUan 1989) in the Cowhng approximation. 

5.1. Outline of Pulsation Calculations 

One possible way for computing nonradial oscillations modes in evolutionary models 
is to first compute such structures with a evolutionary code, and then to choose a subset 
of these models (usually belonging to a predetermined interval in Teff, called hereafter 
Te//-strip) for pulsation analysis. This procedure has been employed in most of the studies 
of adiabatic pulsations in WDs to date (see, e.g., TFW, Bradley, Winget & Wood 1993, 
Bradley 1996, Brassard et al. 1992a, 1992b). 
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Here, we present an alternative way for computing nonradial modes in evolutionary 
sequences. The basic idea is very direct: if the time step in the evolutionary sequence 
is short enough, eigenvalues and eigenfunctions corresponding to a given model should 
strongly resemble those corresponding to the previous one. Then, if we couple the pulsation 
code to the evolutionary code, it is feasible to avoid the scan in u"^ for each model: the 
search of modes is carried out only for the first model inside the relevant Te/j-strip. The 
modes are computed and stored and serve as a trial solution for the next model and so on. 
In this way it is possible economize CPU time, and (more importantly) it is possible to 
follow the changes in the pulsational spectrum due to the evolution of the stellar structure 
in a continuous fashion. 

5.2. Details of Calculations 

Let us briefly describe how our pulsation and evolutionary codes work together. Firstly, 
Te//-strip is chosen, as well as the frequency window to be scanned. The evolutionary 
code computes the model up to the moment when the model reaches the hot edge of 
the Tg/j-strip. Then, the program calls the pulsation routine beginning the scan for 
eigenfrequencies, as described in When a mode is found, the code generates an 
approximate solution for z = 1, ■ ■ ■ , 4 and tu^, which is improved iteratively. Then, such 
solution is tested and, if necessary, meshpoints are redistributed as outlined in and 
iterated to convergence. Before improving the solution, the eigenmodes have been stored on 
the original grid in a common block, with the aim of being employed later as an approximate 
solution for the next stellar model of the sequence. When the computation of the mode has 
been finished, the code begins to look for eigenmodes again, and the process is repeated 
until the given frequency window is fully covered. Thus, we have finished the computation 
of all eigenmodes of the first model belonging to the Te//-strip and the structure of each 
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one is now in the computer memory. Then, the evolutionary code generates the next 
stellar model and the code calls pulsation routine again, but in this occasion the search 
for modes is skipped. Instead of this, stored eigenmodes are taken as input to the general 
Newton - Raphson scheme to approximate the modes of this subsequent stellar model. 
Then, the solution is iterated, and so on. The whole procedure is automatically repeated 
for all evolutionary models inside the Te//-strip. When the star gets outside the Te//-strip, 
computation is finished. 

In terms of CPU, a large fraction of the running time is spent in the search for modes, 
but this is executed only once. The process of relaxation and improvement of the solution 
is faster, and the building of each stellar model is almost instantaneous. 

5.3. Helium White Dwarfs Models 

As mentioned before, we have choose models of hehum WDs for checking the 
efficiency our pulsation code in realistic stellar configurations. To be specific, we computed 
nonradial g- modes for a helium WD model with mass M* = O.SMq, in the range of 
8000ir < Teff < 25000ir. The complete sequence comprise 216 models. Convection, 
present in the outermost layers, is treated with the ML3 version of Mixing Length Theory. 

5.4. Results 

Since one of our interests here is to evaluate the asymptotic behaviour of eigenpcriods, 
we have computed nonradial g- modes for i — 1,2 and 3 with radial order from A; = 1 to 
56. The resulting set of modes is adequate for our purpose, since it covers a broad range in 
periods, enough for enabling the comparison with the predictions of the asymptotic theory 
of Tassoul (1980). Here we employ the mechanical external boundary condition given in 
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Unno, et al. (1989) 



yi 



1 + 



—J- _ 4 _ 











-1/2 + 1/3 
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(30) 



which replaces the first of Eqs. (|T8|). 



For each computed mode, the quantities of interest are the period P^, eigenfunctions 
Ui^i = 1, ■ ■ ■ ,4; kinetic energy K.E., and first order rotation splitting coefficient, Ce^k- In 
addition, we compute the variational period, P^, and the weight function, WF, in the form 
given by Kawaler, Hansen & Winget (1985). Also, for each non-perturbed model we obtain 
the asymptotic spacing of periods (for the same i), APa, given by (TFW, Tassoul 1980) 



APa 



(31) 



where Pq is defined as 



1 iV , 1-1 

— ax 
X 



(32) 



As in TFW, for computing Pq we have ignored the presence of an OCZ, integrating 
from the center to the surface of model but setting = in the OCZ where A^^ < 0. In 
this approximation, we have overestimated the value of AP^, but since the thickness of the 
OCZ (in the radial coordinate) is very small, this is a good approximation. 

The results obtained for £ = 1, 2, 3 are qualitatively very similar. Although pulsations 
have not yet been observed in these stars, in other pulsating WDs the modes with i = 1 
dominate the observations, thus we shall concentrate in £ = 1. 

In Fig. ^ we have plotted the eigenfunction yi for the modes gi, • • •, gs with i = 1, 2, 
for a model with T^jf = 11900 K. From Fig. |^ it is evident that yi has large amplitudes in 
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the whole star, especially in the core. This feature is strongly emphasized when we inspect 
modes with increasing k. This is a remarkable difference in the behaviour of yi comparing 
with the case of oscillating DA and DB WDs. In such WDs yi has a lower central amplitude 
(for modes with £ = 2 in a O.GMq DA WD model see, e.g.. Fig. 14a and 16a of Brassard et 
al. 1992b). 

Fig. 1^ displays the normalized weight function (WF) for gi,- ■ ■,g5 {i = 1) of the same 
model. WF provides a measure of the relative contribution from different zones in the 
model to the period formation. The plot indicates that period formation happens mainly 
in external region of star, where WF is large. However there are important contributions 
from central locations (x ^ 0.2), which is somewhat different to the WF corresponding to 
the case of DA and DB WDs, for which WF displays appreciate values only in the envelope 
(in the "normal" modes; see Figs. 15 and 17 of Brassard et al. 1992b). 

Let us discuss our results concerning periods (P^) and kinetic energies (K.E.). The 
periods here computed show the expected trend in g-modes, with values increasing with 
the order of modes k. Since our calculations covers a wide Te/j-strip, it is possible infer the 
changes of periods during WD cooling. 

It is worth mentioning that we have computed the Brunt - Vaisala frequency as in TFW 
(the Ledoux term B is zero for the whole of model, since it is chemically homogeneous). We 
have found that A^^ decreases as the model cools down, a feature that is very pronounced 
in the core. It is due mainly to the increasing degeneracy in that region, decoupling more 
and more pressure from temperature {xt drops). As consequence eigenperiods increase 
monotonically when Tg// drops below ~ 20000 K, as clearly shown in Fig. 

Next, we consider the behaviour of the kinetic energy (K.E.) of the modes. We have 
found that, as expected in a chemically homogeneous star, K.E. is a smooth function of k 
and thus of the period. In Fig. |^ we show log(K.E.) vs. T^/f for the models included in 
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Fig. 1^. For hotter models it is clear that the more energetic modes are those with lower 
order, because they penetrate deep in the star, where density is high. However, when T^ff 
drops below ~ 20000 K the energy of high order modes (which are mainly concentrated in 
the outer envelope of star) strongly increase. The explanation of this effect (see Brassard et 
al. 1992b in context of a DA WD star) resides in the fact that at such T^ff value the OCZ 
of the object suddenly gets thicker (see the dotted hne in Fig. 17 of Althaus & Benvenuto 
1997). The high order modes that oscillates in the envelope feel gradually the presence of 
convection as the star cools down. Since that in OCZs g- modes become evanescent, such 
modes are forced to get larger amplitudes somewhat below the bottom of the OCZ where 
the density is larger. Since K.E. is proportional to the integral of square of displacement, 
weighted by p, these modes oscillate with larger energies. In sharp contrast to the behaviour 
of high k modes, low order modes are rather insensitive to the thickening of the OCZ. 

Now let us discuss the period spacing of the modes. To our knowledge there is no study 
available in the literature devoted to nonradial pulsations of helium WDs. This situation 
does not allow us to perform a direct comparison of our results. Nevertheless we can have 
a good idea related to the reliability of our pulsation code examining the period spacing 
of consecutive modes of the same i. Tassoul (1980) predicted that the asymptotic period 
spacing for g-modes for a chemically homogeneous object in the adiabatic approximations 
is a constant value AP4 (see Eq. |31|). In Fig. ^ we show the forward period spacing 
APfe = Pk+i-Pk vs. radial order k {^ = I) for models with T^// = 13200, 11800 and 9400 K. 
As reference we also include in this figure the corresponding asymptotic spacing predicted 
by Tassoul (1980) treatment. From Fig. |^ it is clear that the trend of the numerical solution 
is the correct one. 

In order to get an estimation of the accuracy of our treatment, we have calculated the 
variational eigenperiods (P^) for each mode in the formalism presented by Kawaler, et al. 
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(1985). The comparison of variational and numerical eigenperiods for the values of C. here 
considered gives a difference between these treatments which is lower than 1% (see Fig. 
1^ for the case of £ = 1). 

By the way, let us comment on that we have also tested our code running it on the 
same models for which oscillation modes were previously computed. Specifically we have 
employed for such test two carbon - oxygen DA WD models of 0.5 Mq and 0.85 whose 
structure was computed with the WDEC evolutionary code and its vibrational properties 
were previously analyzed (Bradley 1996). In this case MujM^ = 10^^; MHe/M^, = 10^^ and 
Teff ^ 12500K. This is a valuable test for our pulsation code because we could confront our 
results with those of other researchers working on detailed models. Also the fact that these 
models were computed with the same evolutionary code is valuable because differences 
in the resulting periods are due only to the pulsation code and not to differences in e.g. 
between the treatment of the equation of state in our evolutionary code and in WDEC. 
Considering a large amount of modes, the differences between the previously computed 
periods and those we have found has been below ^ 0.1%. Thus, we judge that our code 
produce results accurate enough for our purposes. This comparison has been made possible 
by Paul Bradley who provided us with the DA WD models and the periods of many its 
modes. 



6. Summary and Conclusions 

In this paper we have presented a general code intended for calculating linear, adiabatic, 
nonradial eigenf unctions and eigenfrequencies corresponding to spherically symmetric stars 
by means of a finite difference scheme. This is a very simple strategy based on a slight 
modification of the general Newton - Raphson method for computing stellar evolution 
presented by Kippenhahn, Weigert & Hofmeister (1967) because of the different structure 
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of the matrix to be inverted in finding the iterative corrections. 

As a first test to this code we have apphed it to the simplest case: polytropic spheres. 
As there exist very precise tabulations of the eigenfrequencies for such configurations, we 
have been able to perform a quantitative comparison of our results with the best available 
in the literature (Christensen - Dalsgaard & MuUan 1994). We have found that our code 
is able to compute polytropic eigenfrequencies with an accuracy of 10~^ or much better for 
low k values. 

Then, we have coupled this nonradial oscillation code to our evolutionary code and 
applied it to the computation of the pulsation spectrum of pure, low mass helium white 
dwarfs (WDs). Specifically we considered the case of an object with a mass of 0.3 Mq for 
a wide range of effective temperatures (Te//-strip). In making such computation we have 
employed a new method. Taking a short - enough time step in the evolutionary sequence, 
eigenvalues and eigenf unctions corresponding to a given model strongly resemble those 
corresponding to the previous one. Then, we scan for modes in the first model inside the 
Tg/j-strip. Eigenmodes are computed and stored and serve as a trial solution for the next 
model and so on. In this way it is possible to follow the changes in the the pulsational 
spectrum due to the evolution of the stellar structure in a continuous and fairly detailed 
fashion (see Fig. ||). 

Regarding the computation of g-modes for the 0.3 Mq pure helium WD model we have 
presented eigenfunctions and weight function (WF) for low order k modes, and periods, 
period spacing, kinetic energy (K.E.) and variational periods for a wide range of modes. 
These quantities were computed for the case of £ = 1,2 and 3 but in Figs. || - |^ we have 
presented only the case ^ = 1 for brevity. To our notice this is the first effort in studying 
the pulsation properties of helium WDs. 

Our experience with realistic helium WD models indicates that our numerical code for 
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pulsations work nicely and in particular the solutions we have found have an asymptotic 
period spacing in good agreement with the predictions of the analytical theory of Tassoul 
(1980). This, and the fact that comparison with more complex WD models with layers of 
different internal composition show us that our numerical tool for computing eigenmodes 
work adequately, especially regarding the coupling we have performed between it and our 
evolutionary code. We plan in the near future to apply this new code to the study of 
pulsation of WDs in general. 

The authors would like to warmly acknowledge Paul Bradley for his kindness in 
providing us with numerical models of carbon - oxygen DA WDs together with his pulsation 
calculations and also for the time he spent in doing so. This allowed us to perform key a 
test to our code, such that it made possible to us to be confident with the results the code 
produces. We also want to thank our referee D. Winget for suggestions that allowed us to 
improve the original version of this paper. 
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Fig. 1. — The absolute value of relative difference of the square of the eigenfrequency for 
the case of polytropes with n — 1.5, 3 and 4, for £ = 2 (filled circles), £ — 3 (filled squares) 
and i — A (filled triangles) compared with the very accurate calculations of Christensen - 
Dalsgaard and MuUan (1994). The points corresponding to modes of the same degree are 
connected for clarity. Notice that differences between those sets of computations arc larger 
the larger is the order k of the mode (i.e., in the case of strongly oscillating eigenfunctions) . 
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Fig. 2. — A The eigenfunction yi for the modes gi, ■ ■ ■, gs with £ = 1, corresponding to 
the 0.3 Mq pure hehum WD model with Tg// = 11900K. Notice the large amplitude of such 
modes in the stellar core. B Same as A but with ^ = 2. 

Fig. 3. — The normalized weight function WF corresponding to the same modes included in 

Fig. gA. 

Fig. 4. — A Periods of dipolar modes (£ = 1) with radial order from A; = 1 to A; = 56 and 
B the asymptotic period spacing AP4 predicted by Tassoul (1980) theory as a function of 
the effective temperature. Notice that the periods of high order modes have a very similar 
behaviour when compared with AP4 during WD cooling. 

Fig. 5. — Kinetic energy (K.E.) vs. Tg// for the same modes with £ = 1 included in Fig. ^ 
The unit of K.E. is ergs and we have assumed that the amplitude of yi at surface is one. For 
more details, see text. 

Fig. 6. — Forward period spacing (AP^) vs. radial order, for modes with i = 1 for three 
models with different Tg/j values. Symbols corresponding to modes of the same T^ff are 
connected for clarity. Horizontal lines correspond to the value of asymptotic period spacing 
in each effective temperature according to Tassoul (1980). 

Fig. 7. — The absolute value of relative difference between numerically computed and 
variational eigenperiods for modes with radial order from /c = 1 to A; = 56 and £ = 1 in 
all models analyzed. 
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